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Abstract 

We consider single-file diffusion in an open system with two species A, B of particles. At the 
boundaries we assume different reservoir densities which drive the system into a non-equilibrium 
steady state. As a model we use an one-dimensional two-component simple symmetric exclusion 
process with two different hopping rates Da,Db and open boundaries. For investigating the dy- 
namics in the hydrodynamic limit we derive a system of coupled non-linear diffusion equations for 
the coarse-grained particle densities. The relaxation of the initial density profile is analyzed by 
numerical integration. Exact analytical expressions are obtained for the self-diffusion coefficients, 
which turns out to be length-dependent, and for the stationary solution. In the steady state we find 
a discontinuous boundary-induced phase transition as the total exterior density gradient between 
the system boundaries is varied. At one boundary a boundary layer develops inside which the 
current flows against the local density gradient. Generically the width of the boundary layer and 
the bulk density profiles do not depend on the two hopping rates. At the phase transition line, 
however, the individual density profiles depend strongly on the ratio Da/Db- Dynamic Monte 
Carlo simulation confirm our theoretical predictions. 
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I. INTRODUCTION 



Single-file diffusion occurs when a random motion of particles is confined to a narrow 
channel where they cannot pass each other due to hard-core repulsion. Examples for this 
behaviour include molecules diffusing in the channels of certain zeolites , transport 

of colloidal spheres in one-dimensional channels [^L or the dynamics of "defects" inside the 
"tube" to which entangled polymers are confined 5||. Also biological motors moving along 
macromolecules or microtubuli exhibit single-file motion, see e.g. 0, Q, and references 
therein. Single-file diffusion is an effectively one-dimensional many-body random process 
which exhibits intriguing correlation effects. It is known for a long time that the mean-square 
displacement of a tagged particle grows subdiffusively for late times with the square root of 



time yi, see also 



10| for an insightful theoretical derivation. Also macroscopic stationary 



properties in open systems, driven by an external gradient of the chemical potential shows 
interesting long-range correlations even if the interaction between particles is strictly local 



11 



12, 



13l |. This implies failure of Onsager theory and indicates strong non-equilibrium 



behaviour. 

So far most theoretical treatments of single-diffusion have focussed on indistinguishable 
particles. The reference model for this process is the one-dimensional symmetric exclu- 
sion processes (SEP), a lattice model where hard-core particles hop randomly to nearest- 
neighbour sites, provided the target is empty However, already an experiment aimed 
at verifying the subdiffusive behaviour of a single tracer particle requires to tag one or more 
particles to make them distinguishable. This is done by labelling some particles without 
changing their diffusion properties and thus corresponds to a two-species SEP with identical 
hopping rates. A rigorous treatment of the macroscopic behaviour of such a two-component 
SEP with many tagged particles in terms of coupled non-linear diffusion equations for an 
infinite system was given by Quastel 1^. A finite system with open boundaries was treated 



only very recently 



17|. 



More generally, however, two physically distinct species of particles may simultaneously 
move inside the same channel. In recent measurements a mixture of toluene and propane 
was adsorbed into different zeolites jl^. The authors measured the temperature dependent 
outfiow and noticed a trapping effect due to single-file diffusion: In zeolites with long in- 
tracrystalline channels the stronger adsorbed toluene molecules control due to blocking the 
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outflow of propane. In terms of the SEP the different physical properties of toluene and 
propane require two different hopping rates. Also studying the drift velocity of an entangled 
polymer in gel electrophoresis in the framework of reptation theory necessitate a description 
of the defect dynamics in terms of two distinct species of particles 



Somewhat surprisingly the hydrodynamic theory of two-component lattice gas models 
seems quite undeveloped even though there is some recent progress [21I, Q, Q- In 3 Keil 
et al. review the Maxwell-Stefan theory describing the diffusive behaviour of a binary fluid 
mixture where the total current, i.e. the sum of both species, is zero. The particle-particle 
interaction is taken into account by including a friction between the species being propor- 
tional to the differences in the velocities. However, this approach is purely phenomenological 
and does not provide a link between microscopic properties of the system and the transport 
coeffficients entering the macroscopic equation. In particular, this approach is not adequate 
for single-file diffusion in a finite system where the self-diffusivity of single particles has 
turned out to play an important role |17|. 

In this paper we go beyond Ref. [l^ by considering a two-component symmetric exclusion 
process with different hopping rates (defined in Sec. 2) and by studying not only stationary 
properties but also the time evolution of the particle density and relaxation towards the 
stationary state. To this end, generalizing the arguments put forward in [l^ for identical 
particles, we consider open boundaries where both ends of the chain are connected to particle 
reservoirs with different chemical potentials. We derive in Sec. 3 a set of coupled nonlinear 
diffusion equations and we analyze the time-dependent solutions of these PDE's by numerical 
integration. This is compared with simulation data obtained from dynamical Monte-Carlo 
simulations (DMCS) of the process. In Sec. 4 the stationary solution is obtained analytically 
in closed form and discussed in some detail. Some conclusions are drawn in Sec. 5. The 
exact computation of the self-diffusion coefficients for the two particle species is presented 
in the first appendix. The second appendix provides for self-containedness some technical 
details of that derivation. 



II. TWO-COMPONENT SYMMETRIC EXCLUSION PROCESS 



Following 



171 we consider a one-dimensional lattice with L lattice sites (Fig. d]). Each 



site i can be empty or occupied by a particle of type A or B. Due to hard-core interaction 
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tA/B 12 3 L Ha/B 

FIG. 1: Two-component symmetric exclusion model with open boundaries. Each species of 
particles has its own hopping rates Da,Db- At the boundaries particles are extracted and 

injected with rates as indicated. 

any site carries at most one particle. Particles can hop to nearest neighbour sites (provided 
the target site is empty) with hopping rates Da/b which are the "bare" diffusion coefficients, 
i.e., the proportionally factor in the mean-square displacement if no other particles were 
present. Jump events occur after an exponentially random time which in dynamic Monte 
Carlo simulations (DMCS) is modelled by random sequential update |17|. 

Let ai (bi) count the A (B) particles on site i. Then the densities are the ensemble- 
averaged expectation values of the respective particle occupation numbers: (oj) = PAi'i), 
{bi) = psii). The probability of finding no particle at site i is (fj) = (1 — — hi). When 
we consider a chain with open boundary conditions, particles are injected and removed 
according to the boundary rates oia/b^Ia/b, Pa/b and 5a/b- So the attempt rate for an 
A-particle to enter the system at the left boundary is a a- It leaves the channel at the 
left boundary with rate 7^ as illustrated in Fig. [TJ The other boundary rates are defined 
similarly. 

n 

In terms of a quantum Hamiltonian formalism [15| the probability distribution of the 
system is represented by the probability vector \P{t) > which evolves in time according to 
the master equation 



-H\P(t) > 

where the generator H has the structure 

L-l 

H = bi + bL + ^ hi^i+i- 



(1) 



(2) 



4 = 1 



For an explicit representation of the generator we assign to the three allowed single-site 
states the three basis vectors 



\A >= 11 >= 







|0>= |0>: 



1 

VV 



\B >E 



1 >= 







(3) 
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corresponding to having an A, no particle or B at site i, respectively. A configuration of the 
full lattice is then represented by a tensor product of these basis vectors. To construct H, 
let E^'^ be the 3x3 matrix with one element located at column x and row y equal to one. 
All other elements are zero: (i?^'^)^;, = S^^a^y^b- The operator for annihilation (creation) 
of an A-particle at site k is then = eI''^ (a^ = eI'^) and for annihilation (creation) 
of B is = E^"^ (pl = EI'^). Finally, the number operators are = bk = E^'^, 

Vk = 1 — ttk — bk. In this representation the boundary matrices are 

bi = aA{vi - al) + asivi - bf) + 7A(ai - a^) + 7i?(&i - K) (4) 
bi = Sa{vl - a'l) + 6b{vl - b^) + (3A{aL - a^) + Psibi - b^). (5) 

Hopping in the bulk between site i and z + 1 is generated by 

hi^i+i = DA{aiVi+i + Vitti+i - a^ajj^^ - a+ar^J 

+ DBihvi+i + Vibi+i - bfb-_^-^ - b'bt^-^). (6) 

The model is now well defined. We refer to appendix A for more details on the formalism. 



In ITfl it was shown that with the choice ^ = |^ = /i^ and ^ = |^ = /i^ the process 
has an uncorrelated stationary equilibrium distribution with constant local particle densities 

n = 7A ^ /3a 

F^iiaA_|_aBi|5A_|_iB 

7A ^ 7s /3a /3s 

„ _ 7s _ /o^ 

1 I OA _|_ OS 1 I 5^ I £b ■ 

7A ^ 7s /3a /3s 

Besides giving the equilibrium densities, these two equations above provide a recipe of how 
to translate the picture of inserting and deleting particles on boundary sites into a picture 
of constant reservoirs at the ends, see Fig. [2l One may regard the middle expression in ((71) 
as a left reservoir density of A-particles, located at a virtual site 2 = 0, with which the 
system exchanges particles by jump events with rates Da,b- Similarly, the middle expression 
in ([8]) is a regarded as left reservoir density of -B-particles. At the right boundary one 
adds a virtual reservoir (located at site i = L + 1) with boundary densities p\, p^ given by 
the rightmost expressions in (I7|), (IH|). 

In this picture, jumping from a reservoir site onto the boundary of the system occurs 
proportional to the respective hopping rate and proportional to the single species boundary 



FIG. 2: Three-state symmetric exclusion model with open boundaries - reservoir picture, 
density. This parameterisation satisfies ([7]) and (JHj) if 

aA/B = Da/bPa/b^ 7A/B = DA/Bil ~ Pa- Pb) (9) 

Sa/b = Da/bp\/b^ Pa/b = Da/b{1 -Pa- Pb)- (10) 

The second equality in (JT]), ((HI) is then the equilibrium condition p^^ = p\q of equal 
reservoir densities or chemical potentials Pab — Pab- ^ non-equilibrium setting the 
left reservoir densities p^, are not equal to those on the right boundary {p\, p^). In 
this boundary-driven case the system evolves towards a stationary state with non-vanishing 
particle currents. 

Notice that for periodic boundary conditions the system is highly non-ergodic not only 
because of particle number conservation, but also because of the single-file constraint which 
conserves the sequence of particles. The open system, on the other hand, is ergodic because 
of the violation of the conservation of particles at the boundaries. 



III. RELAXATION 



A. Hydrodynamic limit 

The average densities (oj) and satisfy the equations of motion ^ (aj) = — (ajiJ), 
^ {bi) = — {biH). This provides the Master equations for a single site, 

^ {tti) = Da i{ai-iVi) + {tti+iVi) - (aifj+i) - (aifj.i)) (11) 
^ (bi) = Db {{bi-iVi) + {bi+iVi) - {biVi+i) - {kvi^i)) . (12) 

At this stage we assume an infinite system and do not care about boundary sites. Neverthe- 
less, in this form the equations of motion are not integrable. Replacing the joint probabilities 
by products of expectation values - according to a simple mean field ansatz which has been 
proven useful in other exclusion processes - fails to provide an even qualitatively correct 
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description. However, an exact equation containing no correlators can be obtained from a 
weighted sum 

d [{a,) + = ((^^_^) + (^^^^) _ 2 (a,) + {h,_,) + - 2 (6,)) . (13) 



dt \ Da Db, 
of the individual particle densities. 

The right-hand side contains a second-order difference for both species individually. 
Coarse-graining, i.e., replacing i by the continuous variable x = ^ and introducing pA{x,t)^ 
PB{x,t) as the A, S-particle density at x, transforms (fT3ll (for vanishing lattice constant a) 
into the macroscopic equation 

+ -5.W,t)+P.(x,t)). (14) 



Da Db 

Here we used diffusive rescaling of the time-coordinate. Notice for the hydrodynamic limit 
(fT4ll to be valid we implicitly use good local ergodic properties in the identification of the 
lattice expectation with the coarse-grained deterministc particle density. The weighted den- 
sity on the l.h.s. of (fT4|l and the total density on the r.h.s. reappear in the computations 
below and for notational simplicity we introduce 

p{x, t) = pa{x, t) + pb{x, t)a{x, t) = — + — (15) 

which reduces to a = p/ D for equal single-particle diffusion coefficients Da = Db = D. 

The linear equation (fTSll does not provide a full description of the coarse-grained dynam- 
ics. To achieve this goal we generalize the approach of [l^ and make an ansatz for the 
dynamics of a single particle (of species A or B) localized at position x. This test particle 
acts as a tracer particle in the background of other particles, with a diffusive motion par- 
tially determined by its self-diffusion coefficient. Additionally, the test particle is subject to 
a background drift b caused by the collective evolution of the entire system towards station- 
arity. A priori the self-diffusion coefficient for the two species of particles could be different. 
However, due to the single-file condition the long-time dynamics of a tracer particle is en- 
tirely determined by the background of other particles which is the same for both species. 
Hence we expect identical self-diffusion coefficients Dg This physical picture leads us 



to the ansatz 



dtPA{x, t) = dlDsPA{x, t) - dj){x, t)pA{x, t) (16) 
dtPB{x,t) = dlDspB{x,t) - dxb{x,t)pB{x,t). (17) 
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for the time evolution of the coarse grained particle densities. The drift term b can then be 
determined by using (fT4|) 



b = -d^{Ds(T - p). 
a 



(18) 



In an infinite system the self-diffusion coefficient vanishes, as is indicated by the subdiffu- 
sive nature of single-file diffusion which was proved rigorously for tracer diffusion in the SEP 
in 25 1 . However, as argued in we expect that in a finite system with open boundaries 
correction terms of leading order 1/ L appear. This is confirmed by the exact result 

1 1-p 



a 



(19) 



derived in Appendix [A] for a periodic system. Introducing the two-component local order 
parameter 



(20) 




P = 

finally leads to the coupled nonlinear diffusion equation 

dtp = (dO^p] 

with the diffusion matrix 



(21) 



1 




' Pa 


Pa 1 






1 


a 




{PB 


PB J 





PB 


PA 


Dg 


Db 


PB 


PA 


Da 


Da 



(22) 



This equation provides a complete description of the density dynamics in the hydrodynamic 
limit. For given initial density profiles the solution of (fT6l) . ffTTl) is uniquely deternined. 
Notice that for an infinite system where Dg = Q the diffusion matrix D has a vanishing 
eigenvalue. For Da = Db (|T6l) and (fTTI) reduce to the coupled nonlinear diffusion equations 
proved rigorously in [l^. 



B. Time evolution of density profiles 

We have no analytical expression for the solution of the system (l2Ti l. but numerical 
integration can be performed to great accuracy using standard routines. We first study 
the case of equal hopping rates Da = Db- This case is similar to the tracer diffusion 

8 



profile D^sDg 



profile D;\=Db 



0.6 

0.5 




150 



200 



profile Da=Db 




0.6 
0.5 
0.4 
0.3 
0.2 
0.1 


0.6 
0.5 



b) 






- 








\ 






t=5000 






\\ Pb 




J 

/J 



50 



100 
lattice site 

profile D;\=Db 



150 



200 



200 



0.2 
0.1 





d) 








.^-"^ 


Pa 










t=400000 






Pb 







50 



100 150 
lattice site 



200 



FIG. 3: Density profiles of the y4(i?)-component (upper/lower curve) at diff'erent MC times 
for a system {L = 200) with equal diffusion coefficients. Snapshots are taken at a) t = 0, b) 
t = 5000, c) t = 50000, d) t = 400000. The system relaxes towards the equilibrium state 
given by the reservoir densities p^j^ = 0.3. Symbols denote densities obtained from 
DMCS. Full lines come from numerical integration of the diffusion equation. 



investigated in j27||. Fig. [3] shows the penetration of tracer particles into the system with 
lattice size L = 200 at four different times. Full lines are solutions of (l2Tll obtained from 
numerical integration. Circles indicate the numerical DMCS densities. Initially, the lattice 
is prepared with a uniform distribution of 5-particles and very few A-particles. The smooth 
drop of the initial A-particle profile at the edges was chosen for faster numerical integration. 
The boundary reservoir densities are Pa=Pa — Pb — Pb — For all instances of 

time the numerical integration shows a rather satisfying collapse with the profiles obtained 
from Monte Carlo simulation. Any DMCS in this work was averaged over at least 10000 
realizations. 

It is remarkable that the PDEs describe the intermediate fiattening at the boundary (see 
e.g. Fig. [3b) of the A-particle profile which was observed for the first time in [27||. It is also 
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interesting that the spatial mean density of -B-particles first grows before it relaxes to its 
equilibrium value 0.3. The explanation comes from the hydrodynamic equation ([T4]) which 
states that the total density profile, i.e. the sum of both species, evolves according to an 
ordinary diffusion equation. Hence, the relaxation to equilibrium is proportional to L^, in 
contrast to the much slower relaxation time of order of the individual components. The 
5-particles compensate the A-particle profile until both reach equilibrium. 

As a second example we consider the relaxation of a mixed system with heavy B particles 
and light A-particles where Da/ Db = 5. For Fig. [4] a)-c) we have chosen a model with 
L = 250, and reservoir densities Pa = Pa — Pb — Pb — 0-3 as before. The agreement 
between simulation and integration in the vicinity of the boundary is not as good as for tracer 
particles. Nevertheless, comparing the simulation results for the smaller lattice (L = 100) 
shown in d) with the results for the larger lattice shown in c) suggests that the deviations 
are finite-size effects. We do not have a theory of finite-size effects beyond the leading 1/L- 
correction to the self-diffusion coefficient. Therefore we do not pursue this issue further. 



C. Small diffusion rate for one of the species 

We discuss the case where one of the species, say B, is very heavy and has a vanishing 
hopping rate Db ^ 0. In this case even the open system is not ergodic and the stationary 
state to which the system relaxes depends on the initial distribution of A and 5-particles. B- 
particles just remain at their initial position and the equilibrium A-particle density between 
two 5's is just their local initial density. Let us now assume an initial state without any B- 
particles in the channel and let the reservoir densities be non-vanishing with Pa = Pa ~ Pa 
and pb = Pb = Pb- For the case of strictly vanishing hopping rate Db = both injection 
rates of 5-particles at the boundaries vanish and therefore no 5-particles ever enter the 
system. Then the bulk dynamics follow the rules of the single-species SEP, but an anomaly 
comes from the fact that injection of A-particles is proportional to the A-particle density 
whereas the outgoing rate occurs according to the hole density which involves the reservoir 
density pb 7^ 0. Hence, A-particles do not equilibrate towards pa (as expected from the 
single-species SEP) but reach a higher value: 

Db=o ^ o^A ^ Sa ^ Pa ^23) 
^ aA + lA Sa + f^A I- Pb' 
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FIG. 4: Snapshots of the density profiles for a system with L = 250, Da/ Db = 5 and open 
boundary conditions with reservoir densities p^/g = 0.3. MC times are a) t = 0, b) 
t = 10000, c) t = 50000. Circles denote densities obtained from DMCS. Full lines come 
from numerical integration of the hydrodynamic equation. The stronger discrepancy 
between simulation data and integration shown in Fig. d) for L = 100 and t = 5000 
demonstrates significant finite size effects. 

Notice that any ratio Pa^~^/pa > 1 can be achieved provided ps and Da are large enough. 

This leads to an interesting conclusion for small but non- vanishing < Db Da/L"^- 
Then channel is in a quasi-equilibrium with A-particles (with density fl23l) ) before eventually 
a i?-particle enters the system. Only very slowly both species start to relax to their com- 
mon equilibrium state and one expects the A profile to initially exceed its true equilibrium 
value Pa, up to some maximal value which is bounded by the maximal density f l23l ). This 
overshooting is indeed observed in simulation, see Fig. [5]for the time evolution of the spatial 
mean densities. Initially the system (L = 50) is almost empty and the particles are injected 
and removed at the boundaries according to the reservoir densities p^/^ = 0.3. The ratio of 
the diffusion coefficient is Da/ Db = 30. One can clearly see the excess of A-particles which 
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spatial mean density (D^/Db=30) 
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FIG. 5: Mean density of A and -B-particles as a function of MC steps. The relaxation 
towards the equilibrium values p^/^ = 0.3 was started from an almost empty lattice 

{L = 50). 

disappears only very slowly. 



IV. STEADY STATE BEHAVIOUR 



A. Finite L 



In the steady state the time derivative in the diffusion equation (1211 1 vanishes and from 
dm) we obtain the integration constant 

d,p = -c. (24) 

In terms of the individual stationary particle currents we have c = Ja/Da + Jb/Db- Inte- 
gration yields the stationary total density profile p = Pa + Pb 

(25) 



/ \ - P P 

p[x) = p H ; X 



L 

in terms of the boundary values p^ . In terms of the boundary densities the integration 
constant is therefore given by 

c = . (26) 
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In order to compute the density profile of the A-species we introduce 

Pa{x) DbPa{x) 



h{x) = 



a{x) 



1 _ £b. 

^ Da 



Pa{x) 



(27) 



Given h both pA and ps can be computed from (i25ll . Using the self-diffusion coefficient (1191) . 
the stationarity condition for pA becomes 

^d, [(1 - p)h] + (l + ^) hd^P = -3 A (28) 

with an integration constant ja proportional to the current of A-particles. Using the linearity 
of the total density profile this yields 



^^^^/l'(x) - Ch{x) + J A = 0. 



where the prime denotes the derivative w.r.t. x. 

This differential equation is straightforwardly integrated and one finds 

L 



(29) 



h{x) 



h- - 



1-1 



p — p X 
1- p- L 



(30) 



Here we have introduced the boundary value 



Pa I Pb 
Da ^ Db 



(31) 



Solving for the density profile of A-particle then yields 



X 



Pa{x) = p + (p+ -p )- 



h{x) 



Db+[1- ^ j Hx) 
The integration constant Ja is given in terms of boundary densities by 

. h+il-p-)'-h-{l~p+fp--p^ 
{l-p-f-{l-p^t L 



(32) 



(33) 



Generically this quantity is of order l/L, up to corrections which are exponentially small in 
system size as long as 7^ p+. 

For vanishing reservoir gradient p^ = p^ = p we have c = 0. This yields a linear function 



X 



h{x) = h- - {h- - h+)- 



(34) 
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and nonlinear density profile 

Pa{x) = 
The current of A-particles 



Pa 


1- 




DB^ 
Da J 


Pa 
P . 


+ {p\ 


-Pa) I 




Da J 






Pa- Pa X 
P L 



(35) 



(36) 



is only of order as opposed to the generic 1/L dependence. For the total particle 

current we find 

p(l-p) {p\-Pa) {Da'-D^') _ P, 



J =JA+JB 



(37) 



L2 a+a- 

which is also of order IjL?. We remark that this current is proportional to the boundary 
gradient of the self-diffusivity. 



B. Phase transition in the thermodynamic Hmit 

The expression for the A-particle density derived for finite L is cumbersome. Rather 
interesting behaviour emerges from an analysis of the thermodynamic limit L ^ oo. To 
leading order in 1/L we find 

^ 4- 

for < p+ 



Pi 


P 


--P+ 


Pa I Pb 

Oa^Db 




L 


Pa 


P~ 


--P+ 


Pa I Pb 

Da~^Db 




L 



jA = { „_ „^ (38) 

for p > p^ 

Hence, as one expects, the A-current is always opposite the total reservoir gradient Ap = 
— , but it changes in a non-analytic fashion at = p"^ where it vanishes to leading 

order in 1/L, see (iSGl ). For positive reservoir gradient Ap the amplitude of A-current is 

governed by the reservoir density at right boundary, otherwise by the left boundary. 

A nonanalyticity at Ap = appears also in the behaviour of the density. Using 

lim2,^oo(l — ax/ L)^ = e~"^ yields the following behaviour. 

Case 1, p^ < p^ : Straightforward algebra turns fl30l) into 

h{x) = h+ + {h- - h+)e-^/^ (39) 
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with localization length 



In 



This gives 



PAix) = p + (p+ - p ) 



X 



p - 



{h- 



-xii 



Da 



[h+ + {h- - /i+)e-^'/C] 



(40) 



(41) 



For large distance x ^ ^ from the left boundary the function h approaches a constant 
and the density profile of A-particles becomes proportional to the total density 



Pa{x) 



(42) 



Therefore at large distance from the boundary the density of each species becomes a linear 
function. The exponentially decaying deviation from linearity marks the occurrence of a 
boundary layer of width ^ at the left boundary. Boundary layers in two-component systems 
with open boundaries were previously observed in 
asymptotic space-averaged mean density 
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2611 and in earlier work 



The 



PA 



lim — / dxpA^x) 



takes the value 



PA 



Pa P + P~ 



(43) 



(44) 



p+ 2 

which is independent of Db/Da- 

Case 2, > p^ : The computation for negative reservoir gradient Ap < is analogous and 
all results can be obtained from Case 1 by interchanging (+, — ) and (x, L — x). One finds 

h{x) =h- - {h- - /i+)e-(^~^)/« (45) 

with localization length 



In 



P -P 
1-P+ 



n -1 



(46) 



This corresponds to a boundary layer at the right boundary of the system. For large distance 
L — x^^ from the right boundary the function h approaches a constant and the the density 
profile of A-particles becomes proportional to the total density 

Pa 



pAix) 



P" 



^p(x) 



(47) 
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Correspondingly one has 



p, = ^i^ (48) 



Therefore the mean density as a function of the reservoir gradient Ap has a jump disconti- 
nuity at Ap = 0. In each of the two phases it is controlled by the boundary which does not 
exhibit the boundary layer. 

V. CONCLUSIONS 

The paper is dedicated to a detailed quantitative investigation single-file diffusion with 
two species of particles in an open system. Adapting ideas from the probabilistic hydro- 
dynamic approach and using information from the microscopic master equation lead us to 
derive a nonlinear diffusion equations for the evolution of the macroscopic particle densities. 
We compute the exact self-diffusion coefficients Dg for the two species which are equal and 
of order 1/ L. Hence they vanish in the hydrodynamic limit, in agreement with the subdiffu- 
sive behaviour of single-file diffusion. However, setting = renders the diffusion matrix 
singular (with a vanishing eigenvalue) and the boundary-value problem becomes overdeter- 
mined. Therefore we conclude that Dg is required for regularization. By use of numerical 
integration and Monte Carlo simulation the relaxation for the two species towards equi- 
librium is analysed. The very slow relaxation of the individual particle species is affected 
by a fast relaxation of the sum of both species. The density of the fast species reaches a 
maximum before it starts to relax slowly to its actual equilibrium value. From the diffusive 
nature of the collective behaviour we conclude that the collective relaxation time is of order 

. The 1/L-behaviour of the self-diffusion coefficient then implies single-species relaxation 
times of order . 

The diffusion equations are solved analytically for the stationary case. There is a 
boundary-induced non-equilibrium phase transition of first order. This transition resem- 



bles a similar transition observed in previous work [l7| in so far as (i) this transition occurs 
when there is no gradient in the reservoir densities p^, p~ and (ii) the transition is accompa- 
nied by a discontinuous change of the position of the boundary layer from one boundary to 
the other. However, a detailed analysis of our results reveals some interesting new features. 
Neither the width of the boundary layer nor the stationary bulk density profiles (l42l) depend 
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on the ratio Db/Da inside the two phases. At the phase transition hne Ap = 0, however, 
the individual density profiles are strongly sensitive to the ratio Db/Da- Unlike in the case 
of equal hopping rates the total particle current j does not vanish, even though the total 
density p = Pa + Pb is. constant and equal to the reservoir densities. In terms of currents 
the phase transition occurs when the weighted current c = Ja/Da + Jb/Db vanishes. 

Boundary layers are known to occur also in bulk-driven lattice gases and have been 



3, 



30 



Slfl- They arise as a consequence 



analyzed in some detail in a series of recent papers 
of short-range correlations or because of the occurrence of shocks. They have a microscopic 
width which is finite on lattice scale and hence leads to a boundary discontinuity when 
under hydrodynamic scaling the lattice spacing is sent to zero. This is contrast to the 
boundary layers observed here which remains finite even for vanishing lattice spacing. Since 
there are no shocks in purely boundary-driven system, a possible explanation for the origin 
of the boundary layers could be the presence of long-range correlations in the stationary 
distribution. 
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APPENDIX A: EXACT TWO-SPECIES SELF-DIFFUSION COEFFICIENTS ON 
A FINITE LATTICE 



The self-diffusion coefficient is defined to be proportional to the asymptotic variance of 
the displacement x(t). More precisely, we define 

of a tagged particle in infinite space. For an equilibrium system not driven by external forces 
one has {x) = and the variance reduces to the mean square displacement (x^). For finite 
time one cannot generally expect Ds to be constant (anomalous diffusion), in particular for 



single-file systems with one species of particles Dg vanishes asymptotically as t 2, see 



25|. 



and for a mathematically rigorous proof for the SEP 

For a finite system defined on a ring with periodic boundary conditions the definition 
f lAll ) becomes meaningless. For a lattice model we consider instead the number of jumps x"*" 
in "positive" direction (say, clockwise) and the number of jumps x~ in the opposite direction. 
As displacement we then define the quantity The diffusion coefficient may 

then be defined by (lAlll . In a periodic single-file system with L sites Dg is finite even for 
t — s> 00 and is proportional to 1/L 10|. In the two-component symmetric simple exclusion 
process we one at first sight sight expect two different self-diffusion coefficients Dsa and Dsb 
for the A and B particles, respectively. In this section we show that Dsa = Dsb = Dg by 
an exact calculation for t ^ 00 for a finite system with periodic boundary conditions. The 
calculation is done via a mapping of the exclusion process to the zero range process (ZRP) 
as follows. 

In the lattice gas picture the state of the system with length L is characterized by a set 
of L occupation numbers y = {yi, ...,yL) each of which can take values —1 (B-particle), 1 
(A-particle) or (vancancy). One defines a different lattice by labelling the N particles as 
illustrated in Fig. [6l The individual labels correspond to the site number of the new lattice 
which, therefore, is of length N. Let s = (si, s^v) be the configuration of the new lattice. 
Then Si = 0,1,2,.., counting the number of particles at site i is defined as the number 
of empty sites to the right of particle i in the SEP. This yields a particle system without 
exclusion, named zero range process (ZRP) [3^. A jump of particle i in the SEP to the 
right induces a jump of a particle located on site i of the ZRP to site i — 1. Similarly, a 
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-4-3-2 -10 1 2 3 4 




-4-3-2-1 1 2 3 4 

FIG. 6: Mapping of the SEP with two species of particles to the ZRP. Particles in the SEP 
correspond to sites in the ZRP and interparticle distances correspond to occupation 
numbers. The particle hopping rate in the SEP maps to a hopping rate across bonds as 

indicated. 

jump to the left in the SEP corresponds to a jump from i — 1 to i in the ZRP. The different 
hopping rates Da, Db assigned to a particle in the SEP are translated into different hopping 
rates across bonds in the ZRP. The symmetry in hopping of the SEP is translated into a 
symmetric hopping across a bond. 

In order to proceed with the calculation of the A-particle self-diffusion coefficient we tag an 
A-particle and without loss of generality label it by zero. Let us now track the displacement 
of this particle set x(0) = 0. Then the displacement x{t) for the particle located on the ring 
is the number of jumps to the left minus the number of jumps to the right, performed within 
the time span and t. In terms of the ZRP this number corresponds to the total particle 
current across bond (0, 1). The location of a tagged particle is thus given by counting the 
number of jumps across a fixed bond in the ZRP. 

In order to set up the generator for the ZRP we define the single-site basis vectors for site 
i as follows: An infinite row vector with a '1' at the first entry (zero elsewhere) is assigned 
to having no particle at site i. Provided site i carries one particle, then a vector with '1' 
at the second entry (zero elsewhere) is assigned and so on. The basis vectors for the entire 
chain are as usual composed by a tensor product of the single-site vectors. A hopping event 
from site i to i + 1 is described by the combined action of the creation and annihilation 
matrix and . See appendix B for a representation of the matrices and details of the 
calculation. To each bond {i — 1, i) we need to assign a hopping rate Di = {Da, Db) which 
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corresponds to the hopping rate of the i'th particle in the SEP. This yields the following 
expression for the Hamiltonian: 

HzRP = - ^ A {a^af+i - di + afa^^^ - di+i) (A2) 

i 

The diagonal part of the Hamiltonian which ensures conservation of probability is con- 
structed by the operators di. We remind the reader that the ground state of this Hamilto- 
nian gives the stationary probability distribution which is unique for given total number of 
particles in the zero-range process. 

For purposes that will immediately become clear we also introduce a non-equilibrium 
ZRP where particle hopping across bond (0, 1) is asymmetric, i.e., a particle hops from site 
to site 1 with a rate D^q and with rate D^q^^ it hops from site 1 to site 0. Notice that 
due to periodic boundary conditions site is identical with site A^. We set q = e^/^ where 
E plays the role of a local driving force (measured in units of the thermal factor ksT). The 
generator of this driven process may be written 

HzRp{E) = HzRP + V{E) (A3) 

where 

V{E) = V = -DA{q - l){a^at - do) - DA{q-' - l)«ar - ^^i)- (A4) 

In this system a stationary current j{E) will emerge. In the original SEP this modification 
corresponds to a force acting on the tagged particle. The Einstein relation [35|| then asserts 
that 

Ds=f{0) (A5) 
where the prime denotes the derivative w.r.t. the force E. 



The proof of this relation 



does not seem to be generally known and we outline the 



lere. In order to track the number of jumps over bond (—1,0) we follow the 



main ideas 

strategy of |33| and introduce the counting operators and X", counting clockwise and 
anticlockwise jumps, see also [s^ for a more general description of counting processes. The 
counting operators act locally on an additional subspace with the action 

X_|A; > = 1 > (A6) 
X+|A; > = + 1 > . (A7) 
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Tracking the position of the tagged particle in the SEP is now implemented by the operator 
X with eigenvalue k: X\k >= k\k >. The full Hamiltonian H = Hzrp + V may be split 
into a term Hzrp, and a "perturbation" V with V = —Da{X^ — l)aQ — Da{x^ — l)ao a]"- 
Due to the left-right symmetry of the hopping, {x{t)) = 0. It remains to calculate the 
second moment of the displacement x in the late-time limit 

(x2) = (s|XV^*|P*>. (A8) 

where \P*) is the stationary probability vector. This computation starts by expanding the 
exponential exp{—Ht) into a time-ordered Dyson series with respect to the perturbation V 



^-{HzRP+V)t ^ ^-HzRpt 



1- fdnV{n)+ fdn H dr,v{n)V{T2) 

Jo Jo Jo 



(A9) 



with V{t) = e^^^^'^V^e^^^^^'^. When calculating the n'th moment of x any terms of or- 
der than higher n vanish identically in the matrix element flASj) which follows from the 
commutation relation [X, X^] = ±X'^ and (s|X='= = (s|. 



< >= - < s 



|XV^^«^* [ dTiV{Ti)\P* > 

Jo 

+ < s|a;2e-^^«^* [ dn [ ' dT2V{Ti)V{T2)\P* > (AlO) 
Jo Jo 



and one finds 



= Da \im {{do) + {d,)) - Da j^^ dt {s\{do - di)e-^*(do - di)\P*)^ . (All) 

The proof of the Einstein relation (1A5I) is completed by noting that the same expression in- 
volving the time-integrated current-correlation function appears in the computation of j'(0) 
using standard first-order time-dependent perturbation theory from quantum mechanics. 

For the actual computation of Dg we therefore need j{E) which is computed in the next 
appendix. From (IBIOI) we read off by taking the derivative at E = 

1 

.i=l 

Using M = Na + Nb and (IBllI) this yields in terms of the parameters of the two-component 
SEP 



ea'-ia-'4( 

i=i ^ 

from which (fTOjl follows. 



Pa Pb 
Da Db 



(A13) 
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APPENDIX B: QUANTUM HAMILTONIAN FORMALISM FOR THE ZERO- 
RANGE PROCESS 



The state of a single site i is characterized by the occupation number Sj. Let us first 
confine to the case of non-conserved particle number. For the infinite number of states one 
chooses the vector representation with the single-site basis vectors 

(A /o\ /o^ 



10 >= 






J 



II >= 



1 





12 > 




1 



(Bl) 



The basis vectors representing the state of the entire lattice with M is given by the tensor 
product of the single-site vectors 



\ri >= \so > ®|si > 



\sm-i > 



(B2) 



In this basis the operation of creating (af ) and deleting a particle (a^ ) at site i is repre- 
sented by the matrices 



v 



a, 



J 







1 




V 



(B3) 



The operator ala^_^_^, hence, describes hopping from site i + 1 to i. According to the rules 
described in appendix A, conservation of probability demands to introduce 







v 



7 



(B4) 



The stationary state for Hzrp (1A2P can be obtained by a product ansatz of the form 



IP* >= (l-z) 



(B5) 



22 



where the fugacity z is related to the ZR particle density via z = For space-dependent 
hopping rates the same ansatz works with space-dependent fagacity Zi. This fugacity is 
determined by the stationarity condition 



H\P* >= (B6) 



which together with 



(di) = Zi (B7) 



yields the recursion relation 



J = Di{zi - Zi+,) ty^O (B8) 
= Doiqzo-q-'z^). (B9) 

Here we have assumed periodic boundary conditions with site M = 0. 
This recursion is solved in terms of the free parameter zq by the relation 

j = ^0 u 1 T- (BIO) 



The fugacity Zq fixes the local density at site i = and via (IB7P all other local densities. We 
remark that Zk — Zi = j X]i=i ^ ~ which by the law of large numbers implies a linear 

fugacity profile on macroscopic scale. For q = I (vanishing driving field = 0) we have 
j = and therefore constant fugacities Zk = z = c/ {1 + c). In terms of the particle density 
in the underlying SEP with N = M particles this gives 

z = l-p. (Bll) 
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